
#### Preliminaries ####

R.version$version.string

# clear workspace
rm(list=ls())

setwd("W:/paper zeitschrift für Menschenrechte")

# load packages
library(here)
library(tidyverse)
library(rworldmap) # maps
library(magrittr)
library(reshape2)
library(ggpubr)
library(leaflet)
library(plotly)
library(gapminder)
library(countrycode)
library(ggrepel)
library(ggthemes)
library(modelsummary)
library(rcartocolor)


#---------------------------------------------------#
#--------------------load data ---------------------#
#---------------------------------------------------#

# 1. Load V-Dem data (version 12)

# set you working directory here #

#setwd()

vdem_full_v13 <- readRDS("C:/Users/ba72loko/projects/data/V-Dem 13 Final/V-Dem-CY-FullOthers_R_v13/V-Dem-CY-Full+Others-v13.rds")


vdem_subset <- vdem_full_v13 %>%
  dplyr::select(country_name, year, country_id, starts_with("v2x_polyarchy"),
                starts_with("v2cafres"), starts_with("v2cafexch"), 
                starts_with("v2cainsaut"), starts_with("v2casurv"), 
                starts_with("v2clacfree"))

vdem <- vdem_full_v13

#### load Political terror Scale data ####

pts_data <- haven::read_dta("data/PTS-2022.dta")

pts_data$country_id <- countrycode::countrycode(source= pts_data$COW_Code_A, "cowc", "vdem", warn = TRUE)

pts_data <- pts_data %>%
  rename(year = Year)

#### Load CIRI Human Rights data ####

ciri_data <- load("data/CIRI_data.RData")

ciri_data <- x
rm(x)

ciri_data$country_id <- countrycode::countrycode(source= ciri_data$COW, "cown", "vdem", warn = TRUE)

ciri_data <- ciri_data %>%
  rename(year = YEAR)

#---------------------------------------------------#
#-----Human Rights Data and Democracy --------------#
#---------------------------------------------------#
prism_colors <- carto_pal(name = "Prism")


vdem_subset_figure_1 <- vdem_full_v13 %>%
  filter(year>=1900) %>%
  dplyr::select(country_name, country_id, year, v2x_clphy, v2x_polyarchy, v2x_libdem,v2mecenefm, v2meharjrn) %>%
  mutate(v2x_clphy_invert = scales::rescale(v2x_clphy, to = c(1,0), from = range(c(0.0140, 0.9900), na.rm = T)))

summary(vdem_subset_figure_1$v2x_clphy)
summary(vdem_subset_figure_1$v2x_clphy_invert)



# v2mecenefm: Government censorship effort
# v2x_clphy: Physical Integrity Index
# v2meharjrn: Harassment on journalist 


ggplot(data = vdem_subset_figure_1, aes(x= v2x_polyarchy, y = v2x_clphy_invert)) +
  geom_point(alpha = 0.2) +
  geom_smooth(method = "gam", formula = y ~ s(x, bs = "cs"), color = prism_colors[c(8)], fill = prism_colors[c(8)]) +
  theme_pubr() +
  labs(x = "Electoral Democracy Index", 
       y = "Physcial Integrity Index (reversed)") 

ggsave("outputs/Expertendaten_empirisch_Demokratie.png", dpi = 600, width = 10, height = 10, units = "cm")


ggplot(data = vdem_subset_figure_1, aes(x= v2mecenefm, y = v2x_clphy_invert)) +
  geom_point(alpha = 0.2) +
  geom_smooth(method = "gam", formula = y ~ s(x, bs = "cs"), color = prism_colors[c(8)], fill = prism_colors[c(8)]) +
  theme_pubr() +
  labs(x = "Government Censorship Effort", 
       y = "Physcial Integrity Index (reversed)") 

ggsave("outputs/Expertendaten_empirisch_Government_censorship.png", dpi = 600, width = 10, height = 10, units = "cm")


ggplot(data = vdem_subset_figure_1, aes(x= v2meharjrn, y = v2x_clphy_invert)) +
  geom_point(alpha = 0.2) +
  geom_smooth(method = "gam", formula = y ~ s(x, bs = "cs"), color = prism_colors[c(8)], fill = prism_colors[c(8)]) +
  theme_pubr() +
  labs(x = "Harassment on Journalist ", 
       y = "Physcial Integrity Index (reversed)") 

ggsave("outputs/Expertendaten_empirisch_Harassment_Journalist.png", dpi = 600, width = 10, height = 10, units = "cm")


# Standard Based Data ##

vdem_subset_figure_2 <- vdem_full_v13 %>%
  left_join(pts_data, by = c("country_id", "year")) %>%
  filter(year >=1990) %>%
  dplyr::select(country_name, country_id, year, v2x_clphy, v2x_polyarchy, v2x_libdem,v2mecenefm, v2meharjrn, 
                starts_with("PTS_"), starts_with("NA_S")) %>%
  mutate(v2x_clphy_invert = scales::rescale(v2x_clphy, to = c(1,0), from = range(c(0.0140, 0.9900), na.rm = T)))
summary(vdem_subset_figure_2$PTS_A)

# PTS_A: Amnesty International Data 
# PTS_H: Human Rights Data 
# PTS_S

# v2mecenefm: Government censorship effort
# v2x_clphy: Physical Integrity Index
# v2meharjrn: Harassment on journalist 

#### Human Rights Data ####

ggplot(data = vdem_subset_figure_2, aes(x= v2x_polyarchy, y = PTS_H)) +
  geom_point(alpha = 0.2) +
  geom_smooth(method = "loess", color = prism_colors[c(8)], fill = prism_colors[c(8)]) +
  theme_pubr() +
  labs(x = "Electoral Democracy Index", 
       y = "Political terrors Scale: HRW") 

ggsave("outputs/Ereignisdaten_empirisch_HRW_Demokratie.png", dpi = 600, width = 10, height = 10, units = "cm")


F_A1.1 <- ggplot(data = subset(vdem_subset_figure_2 %>% drop_na(PTS_H)), aes(x= v2mecenefm, y = PTS_H, color = as.factor(PTS_H), fill = as.factor(PTS_H))) +
  geom_point(position = position_jitter(width = .15), size = .25)+
  geom_boxplot(outlier.shape = NA, alpha = 0.3, width = .1, colour = "BLACK") +
  theme_pubr() +
  guides(fill = FALSE, colour = FALSE) +
  labs(x = "Government Censorship Effort", 
       y = "Political terrors Scale: HRW") +
  scale_color_carto_d("prism")+
  scale_fill_carto_d("prism")

ggsave("outputs/Ereignisdaten_empirisch_HRW_Government_censorship.png", dpi = 600, width = 10, height = 10, units = "cm")


ggplot(data = vdem_subset_figure_2, aes(x= v2meharjrn, y = PTS_H)) +
  geom_point(alpha = 0.2) +
  geom_smooth(method = "loess",  color = prism_colors[c(8)], fill = prism_colors[c(8)]) +
  theme_pubr() +
  labs(x = "Harassment on Journalist ", 
       y = "Political terrors Scale: HRW") 

ggsave("outputs/EEreignisdaten_empirisch_HRW_Harassment_Journalist.png", dpi = 600, width = 10, height = 10, units = "cm")

#### Amnesty International Data ####

ggplot(data = vdem_subset_figure_2, aes(x= v2x_polyarchy, y = PTS_A)) +
  geom_point(alpha = 0.2) +
  geom_smooth(method = "loess", color = prism_colors[c(8)], fill = prism_colors[c(8)]) +
  theme_pubr() +
  labs(x = "Electoral Democracy Index", 
       y = "Political terrors Scale: AI") 

ggsave("outputs/Ereignisdaten_empirisch_AI_Demokratie.png", dpi = 600, width = 10, height = 10, units = "cm")


F_A1.2 <-ggplot(data = subset(vdem_subset_figure_2 %>% drop_na(PTS_A)), aes(x= v2mecenefm, y = PTS_A, color = as.factor(PTS_A), fill = as.factor(PTS_A))) +
  geom_point(position = position_jitter(width = .15), size = .25)+
  geom_boxplot(outlier.shape = NA, alpha = 0.3, width = .1, colour = "BLACK") +
  labs(x = "Government Censorship Effort", 
       y = "Political terrors Scale: AI") +
  guides(fill = FALSE, colour = FALSE) +
  theme_pubr() + 
  scale_color_carto_d("prism")+
  scale_fill_carto_d("prism")

ggsave("outputs/Ereignisdaten_empirisch_AI_Government_censorship.png", dpi = 600, width = 10, height = 10, units = "cm")


ggplot(data = vdem_subset_figure_2, aes(x= v2meharjrn, y = PTS_A)) +
  geom_point(alpha = 0.2) +
  geom_smooth(method = "loess",  color = prism_colors[c(8)], fill = prism_colors[c(8)]) +
  theme_pubr() +
  labs(x = "Harassment on Journalist ", 
       y = "Political terrors Scale: AI") 

ggsave("outputs/Ereignisdaten_empirisch_AI_Harassment_Journalist.png", dpi = 600, width = 10, height = 10, units = "cm")


F_A1.3 <- ggplot(data = vdem_subset_figure_2 %>% drop_na(PTS_S), aes(x= v2mecenefm, y = PTS_S, color = as.factor(PTS_S), fill = as.factor(PTS_S))) +
  geom_point(position = position_jitter(width = .15), size = .25)+
  geom_boxplot(outlier.shape = NA, alpha = 0.3, width = .1, colour = "BLACK") +
  labs(x = "Government Censorship Effort", 
       y = "Political terrors Scale: State Department") +
  guides(fill = FALSE, colour = FALSE) +
  theme_pubr() + 
  scale_color_carto_d("prism")+
  scale_fill_carto_d("prism")



##### CIRI Data ####

vdem_subset_figure_3 <- vdem_full_v13 %>%
  left_join(ciri_data, by = c("country_id", "year")) %>%
  filter(year >=1900) %>%
  dplyr::select(country_name, country_id, year, v2x_clphy, v2x_polyarchy, v2x_libdem,v2mecenefm, v2meharjrn, 
               PHYSINT, DISAP, KILL, TORT, ) %>%
  mutate(v2x_clphy_invert = scales::rescale(v2x_clphy, to = c(1,0), from = range(c(0.0140, 0.9900), na.rm = T))) %>%
  drop_na(PHYSINT)

summary(vdem_subset_figure_3$PHYSINT)


F_A1.4 <- ggplot(data = vdem_subset_figure_3, aes(x= v2mecenefm, y = as.factor(PHYSINT), color = as.factor(PHYSINT), fill = as.factor(PHYSINT))) +
  geom_point(position = position_jitter(width = .15), size = .25)+
  geom_boxplot(outlier.shape = NA, alpha = 0.3, width = .1, colour = "BLACK") +
  theme_pubr() +
  labs(x = "Government Censorship Effort", 
       y = "CIRI: Physical Integrity") +
  guides(fill = FALSE, colour = FALSE) +
  theme_pubr() + 
  scale_color_carto_d("prism")+
  scale_fill_carto_d("prism")

library(patchwork)


(F_A1.1 + F_A1.2)/ (F_A1.3 + F_A1.4) +
  scale_color_grey() +
  scale_fill_grey()

ggsave("outputs/SA_A1.png", dpi = 600, width = 15, height = 20, units = "cm")


#---------------------------------------------------#
#----- Figure 3: Human Rights Indicators V-Dem -----#
#---------------------------------------------------#

vdem_subset <- vdem_full_v13 %>%
  dplyr::select(country_name, country_id, year, starts_with("v2clpolcl"), starts_with("v2clacjust"), starts_with("v2x_clphy"), 
                starts_with("v2cltort"), starts_with("v2xcl_dmove"), starts_with("v2xcl_prpty"), starts_with("v2xcl_disc"),
                starts_with("v2x_freexp_altinf"), starts_with("v2pehealth"), starts_with("v2peedueq"), e_wb_pop) %>%
  filter(year >= 1960) %>%
  group_by(country_id) %>%
  fill(country_name, e_wb_pop, .direction = "down")

summary(vdem_subset$e_wb_pop)

# v2clpolcl: Political group equality in respect for civil liberties
# v2clacjust: Social class equality in respect for civil liberty
# v2x_clphy: Physical Violence Index
# v2cltort: Freedom from tortur
# v2xcl_dmove: Freedom of domestic movement
# v2xcl_prpty: property rights
# v2xcl_disc: Freedom of Discussion
# v2x_freexp_altinf: Freedom of Expression and Alternative Sources of Information Index
# v2pehealth: Health equality
# v2peedueq: Educational equality


e_wb_pop_year <- vdem_subset %>%
  filter(year >= 1960) %>%
  group_by(year) %>%
  summarize(e_wb_pop_world_sum = sum(e_wb_pop, na.rm = T))

vdem_subset_pop <- vdem_subset %>%
  left_join(e_wb_pop_year, by = c("year")) %>%
  mutate(pop_weight = e_wb_pop / e_wb_pop_world_sum) %>%
  drop_na(pop_weight) %>%
  group_by(year) %>%
  summarize(v2x_clphy_mean =  weighted.mean(v2x_clphy, pop_weight, na.rm = T), 
            v2x_clphy_codehigh =  weighted.mean(v2x_clphy_codehigh, pop_weight, na.rm = T),
            v2x_clphy_codelow =  weighted.mean(v2x_clphy_codelow, pop_weight, na.rm = T), 
            v2xcl_dmove_mean =  weighted.mean(v2xcl_dmove, pop_weight, na.rm = T), 
            v2xcl_dmove_codehigh =  weighted.mean(v2xcl_dmove_codehigh, pop_weight, na.rm = T),
            v2xcl_dmove_codelow =  weighted.mean(v2xcl_dmove_codelow, pop_weight, na.rm = T),
            v2xcl_disc_mean =  weighted.mean(v2xcl_disc, pop_weight, na.rm = T), 
            v2xcl_disc_codehigh =  weighted.mean(v2xcl_disc_codehigh, pop_weight, na.rm = T),
            v2xcl_disc_codelow =  weighted.mean(v2xcl_disc_codelow, pop_weight, na.rm = T), 
            v2pehealth_osp_mean =  weighted.mean(v2pehealth_osp, pop_weight, na.rm = T), 
            v2pehealth_osp_codehigh =  weighted.mean(v2pehealth_osp_codehigh, pop_weight, na.rm = T),
            v2pehealth_osp_codelow =  weighted.mean(v2pehealth_osp_codelow, pop_weight, na.rm = T)) %>%
  mutate(averaged = "by Population")


vdem_subset_country <- vdem_subset %>%
  group_by(year) %>%
  summarize(v2x_clphy_mean =  mean(v2x_clphy, na.rm = T), 
            v2x_clphy_codehigh = mean(v2x_clphy_codehigh, na.rm = T),
            v2x_clphy_codelow = mean(v2x_clphy_codelow, na.rm = T), 
            v2xcl_dmove_mean = mean(v2xcl_dmove, na.rm = T), 
            v2xcl_dmove_codehigh = mean(v2xcl_dmove_codehigh, na.rm = T),
            v2xcl_dmove_codelow = mean(v2xcl_dmove_codelow, na.rm = T),
            v2xcl_disc_mean = mean(v2xcl_disc, na.rm = T), 
            v2xcl_disc_codehigh = mean(v2xcl_disc_codehigh, na.rm = T),
            v2xcl_disc_codelow = mean(v2xcl_disc_codelow, na.rm = T), 
            v2pehealth_osp_mean = mean(v2pehealth_osp, na.rm = T), 
            v2pehealth_osp_codehigh = mean(v2pehealth_osp_codehigh, na.rm = T),
            v2pehealth_osp_codelow = mean(v2pehealth_osp_codelow, na.rm = T)) %>%
  mutate(averaged = "by Country")

vdem_subset_figure <- rbind(vdem_subset_country, vdem_subset_pop)

prism_colors <- carto_pal(name = "Prism")


F3.1 <- ggplot(vdem_subset_figure, aes(x =year, y = v2x_clphy_mean, color = averaged, fill = averaged)) +
  geom_ribbon(aes(ymin=v2x_clphy_codelow, ymax= v2x_clphy_codehigh), alpha = 0.3) +
  geom_line(size =1.05) +
  theme_pubr() +
  ylim(0,1) +
  scale_x_continuous(breaks = seq(1960, 2020, by = 10)) +
  labs(y = "Physical Violence Index",
       color = "", x = "Year",
       size = "",
       fill = "",
       linetype = "", 
       shape = "",
       title = "Art. 3 Recht auf Leben und Freiheit",
       subtitle = "Physical Violence Index",
       caption = "") +
  scale_colour_manual(values = c("#696969", "#D3D3D3")) +
  scale_fill_manual(values = c("#696969", "#D3D3D3"))
  
F3.2 <- ggplot(vdem_subset_figure, aes(x =year, y = v2xcl_dmove_mean, color = averaged, fill = averaged)) +
  geom_ribbon(aes(ymin=v2xcl_dmove_codelow, ymax= v2xcl_dmove_codehigh), alpha = 0.3) +
  geom_line(size =1.05) +
  theme_pubr() +
  ylim(0,1) +
  scale_x_continuous(breaks = seq(1960, 2020, by = 10)) +
  labs(y = "Freedom of domestic movement",
       color = "", x = "Year",
       size = "",
       fill = "",
       linetype = "", 
       shape = "",
       title = "Art. 13 Recht auf Freizügigkeit",
       subtitle = "Freedom of domestic movement",
       caption = "") +
  scale_colour_manual(values = c("#696969", "#D3D3D3")) +
  scale_fill_manual(values = c("#696969", "#D3D3D3"))


F3.3 <- ggplot(vdem_subset_figure, aes(x =year, y = v2xcl_disc_mean, color = averaged, fill = averaged)) +
  geom_ribbon(aes(ymin=v2xcl_disc_codelow, ymax= v2xcl_disc_codehigh), alpha = 0.3) +
  geom_line(size =1.05) +
  theme_pubr() +
  ylim(0,1) +
  scale_x_continuous(breaks = seq(1960, 2020, by = 10)) +
  labs(y = "Freedom of discussion",
       color = "", x = "Year",
       size = "",
       fill = "",
       linetype = "", 
       shape = "",
       title = "Art. 19 Meinungs- und Informationsfreiheit",
       subtitle = "Freedom of discussion",
       caption = "") +
  scale_colour_manual(values = c("#696969", "#D3D3D3")) +
  scale_fill_manual(values = c("#696969", "#D3D3D3"))


F3.4 <- ggplot(vdem_subset_figure, aes(x =year, y = v2pehealth_osp_mean, color = averaged, fill = averaged)) +
  geom_ribbon(aes(ymin=v2pehealth_osp_codelow, ymax= v2pehealth_osp_codehigh), alpha = 0.3) +
  geom_line(size =1.05) +
  scale_x_continuous(breaks = seq(1960, 2020, by = 10)) +
  theme_pubr() +
  ylim(0,4) +
  labs(y = "Health equality",
       color = "", x = "Year",
       size = "",
       fill = "",
       linetype = "", 
       shape = "",
       title = "Art. 25 Recht auf Wohlfahrt",
       subtitle = "Health equality",
       caption = "") +
  scale_colour_manual(values = c("#696969", "#D3D3D3")) +
  scale_fill_manual(values = c("#696969", "#D3D3D3"))


ggarrange(F3.1, F3.2, F3.3, F3.4, ncol=2, nrow=2, common.legend = TRUE, legend="bottom", labels = c("A", "B", "C", "D"))

ggsave("outputs/Figure_3_Main_Paper.png", dpi = 600, width = 25, height = 30, units = "cm")

#### Abb. 4 #####

#### World Average ####

vdem_fig_world <- vdem_full_v13 %>%
  dplyr::select(country_name, year, e_regionpol_6C, v2xca_academ, v2xca_academ_codehigh, v2xca_academ_codelow) %>%
  filter(year >= 1900) %>%
  group_by(year) %>%
  summarize(v2xca_academ_mean = mean(v2xca_academ, na.rm = T), 
            v2xca_academ_codehigh = mean(v2xca_academ_codehigh, na.rm = T),
            v2xca_academ_codelow = mean(v2xca_academ_codelow, na.rm = T)) %>%
  mutate(averaged = "by Country")


vdem_fig_world_pop <- vdem_full_v13 %>%
  group_by(country_id) %>%
  fill(e_wb_pop, .direction = "down") %>%
  ungroup() %>%
  left_join(e_wb_pop_year, by = c("year")) %>%
  mutate(pop_weight = e_wb_pop / e_wb_pop_world_sum) %>%
  drop_na(pop_weight) %>%
  group_by(year) %>%
  summarize(v2xca_academ_mean =  weighted.mean(v2xca_academ, pop_weight, na.rm = T), 
            v2xca_academ_codehigh =  weighted.mean(v2xca_academ_codehigh, pop_weight, na.rm = T),
            v2xca_academ_codelow =  weighted.mean(v2xca_academ_codelow, pop_weight, na.rm = T)) %>%
  mutate(averaged = "by Population")

vdem_subset_figure <- rbind(vdem_fig_world, vdem_fig_world_pop)


p1 <- ggplot(vdem_subset_figure, aes(x =year, y = v2xca_academ_mean, color = averaged, fill = averaged)) +
  geom_ribbon(aes(ymin=v2xca_academ_codelow, ymax= v2xca_academ_codehigh), alpha = 0.3) +
  geom_line(size =1.05) +
  theme_pubr() +
  ylim(0,1) +
  scale_x_continuous(breaks = seq(1900, 2020, by = 10)) +
  labs(y = "Academic Freedom Index",
       color = "", x = "Year",
       size = "",
       fill = "",
       linetype = "", 
       shape = "",
       title = "Wissenschaftsfreiheit: weltweite Entwicklung",
       subtitle = "",
       caption = "") +
  scale_colour_manual(values = c("#696969", "#D3D3D3")) +
  scale_fill_manual(values = c("#696969", "#D3D3D3"))


vdem_subset_sum <- vdem_full_v13 %>%
  group_by(year) %>%
  summarize(number_of_countries = sum(v2cauni, na.rm = T))%>%
  filter(year >=1900)

p2 <- ggplot(vdem_subset_sum, aes(x = year, y= number_of_countries)) +
  geom_bar(stat="identity", position="stack", width = 0.7) +
  theme_pubr() +
  scale_x_continuous(breaks = seq(round(min(vdem_subset_sum$year) / 10) * 10, round(max(vdem_subset_sum$year) / 10) * 10, 10)) +
  labs(y = "Countries",
       color = "", x = "Year",
       size = "",
       fill = "",
       shape = "",
       title = "Anzahl der Länder",
       subtitle = "",
       caption = "") +
  theme_pubr()

## Plot 3 ##

vdem_subset_p3 <- vdem_full_v13 %>%
  dplyr::select(country_text_id, country_name, year, v2xca_academ, v2xca_academ_codelow, v2xca_academ_codehigh,  
                e_regionpol_6C) %>%
  filter(year %in% c(1966, 1992)) %>%
  mutate(e_regionpol_6C = case_when(e_regionpol_6C ==1 ~ "Eastern Europe\nand Central Asia", 
                                    e_regionpol_6C ==2 ~ "Latin America and\nthe Caribbean", 
                                    e_regionpol_6C ==3 ~ "The Middle East and\nNorth Africa", 
                                    e_regionpol_6C ==4 ~ "Sub-Saharan Africa", 
                                    e_regionpol_6C ==5 ~ "Western Europe and\nNorth America", 
                                    e_regionpol_6C ==6 ~ "Asia and Pacific"))

vdem_subset_p3_1992 <- filter(vdem_subset_p3, year == 1992)
vdem_subset_p3_1966 <- filter(vdem_subset_p3, year == 1966)
var <- "v2xca_academ"

vdem_subset_point <- vdem_subset_p3 %>%
  dplyr::select(country_text_id, country_name, year, v2xca_academ,e_regionpol_6C) %>%
  pivot_wider(names_from = year, values_from = v2xca_academ,  names_prefix = "v2xca_academ_") 

vdem_subset_codelow <- vdem_subset_p3 %>%
  dplyr::select(country_text_id, country_name, year, v2xca_academ_codelow,e_regionpol_6C) %>%
  pivot_wider(names_from = year, values_from = v2xca_academ_codelow,  names_prefix = "v2xca_academ_codelow_") %>%
  dplyr::select(-c("country_name", "e_regionpol_6C"))

vdem_subset_codehigh <- vdem_subset_p3 %>%
  dplyr::select(country_text_id, country_name, year, v2xca_academ_codehigh,e_regionpol_6C) %>%
  pivot_wider(names_from = year, values_from = v2xca_academ_codehigh,  names_prefix = "v2xca_academ_codehigh_") %>%
  dplyr::select(-c("country_name", "e_regionpol_6C"))

vdem_subset <- vdem_subset_point %>%
  left_join(vdem_subset_codelow, by = c("country_text_id")) %>%
  left_join(vdem_subset_codehigh, by = c("country_text_id")) %>%
  mutate(diff_1966_1992 = 0) %>% 
  mutate(diff_1966_1992 = ifelse(v2xca_academ_codelow_1992 > v2xca_academ_codehigh_1966 & v2xca_academ_1992 - v2xca_academ_1966> 0.1,
                                 1, 0), 
         diff_1966_1992 = ifelse(v2xca_academ_codehigh_1992 < v2xca_academ_codelow_1966 & v2xca_academ_1992 - v2xca_academ_1966 < -0.1,
                                 1, diff_1966_1992))

table(vdem_subset$diff_1966_1992)

p3 <-  ggplot(vdem_subset, aes(x = v2xca_academ_1966, y = v2xca_academ_1992, color = e_regionpol_6C, shape = e_regionpol_6C)) +
  geom_point(size = 2) +
  geom_abline(intercept = 0, slope = 1, size = 0.5)+
  geom_text_repel(data = subset(vdem_subset, diff_1966_1992 ==1), aes(x = v2xca_academ_1966,
                                                                      y = v2xca_academ_1992, 
                                                                      label = country_name), 
                  show_guide  = FALSE) +
  theme_bw() + 
  scale_colour_manual(values = c("#DCDCDC", "#D3D3D3", "#A9A9A9", "#696969","#708090", "#000000")) +
  scale_fill_manual(values = c("#DCDCDC", "#D3D3D3", "#A9A9A9", "#696969","#708090", "#000000")) +
  scale_shape_manual(values=c(15, 16, 17, 18, 24, 25))+
  labs(y = "Academic Freedom Index 1992",
       x = "Academic Freedom Index 1966",
       color = "",
       size = "",
       label = "",
       shape = "",
       title = "1966-1992 Vergleich",
       subtitle = "",
       caption = "") +
  guides(color = guide_legend(override.aes = list(size = 3))) +
  theme(legend.position="bottom")

## Plot 4 ##

vdem_subset_p4 <- vdem_full_v13 %>%
  dplyr::select(country_text_id, country_name, year, v2xca_academ, v2xca_academ_codelow, v2xca_academ_codehigh,  
                e_regionpol_6C) %>%
  filter(year %in% c(1992, 2022)) %>%
  mutate(e_regionpol_6C = case_when(e_regionpol_6C ==1 ~ "Eastern Europe\nand Central Asia", 
                                    e_regionpol_6C ==2 ~ "Latin America and\nthe Caribbean", 
                                    e_regionpol_6C ==3 ~ "The Middle East and\nNorth Africa", 
                                    e_regionpol_6C ==4 ~ "Sub-Saharan Africa", 
                                    e_regionpol_6C ==5 ~ "Western Europe and\nNorth America", 
                                    e_regionpol_6C ==6 ~ "Asia and Pacific"))

vdem_subset_p4_2022 <- filter(vdem_subset_p4, year == 2022)
vdem_subset_p4_1992 <- filter(vdem_subset_p4, year == 1992)
var <- "v2xca_academ"

vdem_subset_point <- vdem_subset_p4 %>%
  dplyr::select(country_text_id, country_name, year, v2xca_academ,e_regionpol_6C) %>%
  pivot_wider(names_from = year, values_from = v2xca_academ,  names_prefix = "v2xca_academ_") 

vdem_subset_codelow <- vdem_subset_p4 %>%
  dplyr::select(country_text_id, country_name, year, v2xca_academ_codelow,e_regionpol_6C) %>%
  pivot_wider(names_from = year, values_from = v2xca_academ_codelow,  names_prefix = "v2xca_academ_codelow_") %>%
  dplyr::select(-c("country_name", "e_regionpol_6C"))

vdem_subset_codehigh <- vdem_subset_p4 %>%
  dplyr::select(country_text_id, country_name, year, v2xca_academ_codehigh,e_regionpol_6C) %>%
  pivot_wider(names_from = year, values_from = v2xca_academ_codehigh,  names_prefix = "v2xca_academ_codehigh_") %>%
  dplyr::select(-c("country_name", "e_regionpol_6C"))

vdem_subset <- vdem_subset_point %>%
  left_join(vdem_subset_codelow, by = c("country_text_id")) %>%
  left_join(vdem_subset_codehigh, by = c("country_text_id")) %>%
  mutate(diff_1992_2022 = 0) %>% 
  mutate(diff_1992_2022 = ifelse(v2xca_academ_codelow_2022 > v2xca_academ_codehigh_1992 & v2xca_academ_2022 - v2xca_academ_1992> 0.1,
                                 1, 0), 
         diff_1992_2022 = ifelse(v2xca_academ_codehigh_2022 < v2xca_academ_codelow_1992 & v2xca_academ_2022 - v2xca_academ_1992 < -0.1,
                                 1, diff_1992_2022))

table(vdem_subset$diff_1992_2022)

p4 <-  ggplot(vdem_subset, aes(x = v2xca_academ_1992, y = v2xca_academ_2022, color = e_regionpol_6C, shape = e_regionpol_6C)) +
  geom_point(size = 2) +
  geom_abline(intercept = 0, slope = 1, size = 0.5)+
  geom_text_repel(data = subset(vdem_subset, diff_1992_2022 ==1), aes(x = v2xca_academ_1992,
                                                                      y = v2xca_academ_2022, 
                                                                      label = country_name), 
                  show_guide  = FALSE) +
  theme_bw() + 
  scale_colour_manual(values = c("#DCDCDC", "#D3D3D3", "#A9A9A9", "#696969","#708090", "#000000")) +
  scale_shape_manual(values=c(15, 16, 17, 18, 24, 25))+
  labs(y = "Academic Freedom Index 2022",
       x = "Academic Freedom Index 1992",
       color = "",
       size = "",
       label = "",
       shape = "",
       title = "1992-2022 Vergleich",
       subtitle = "",
       caption = "") +
  guides(color = guide_legend(override.aes = list(size = 3))) +
  theme(legend.position="bottom")

#### Combine Plots ####

p3_4 <- ggarrange(p3, p4,align = "h",
                           labels= c("C", "D"), common.legend = T)



Figure_4 <- cowplot::plot_grid(cowplot::plot_grid(p1, p2,  
                                      ncol = 1, align = "v", 
                                      rel_heights = c(1, 0.6), labels = c('A', 'B')), 
                   p3_4, 
                   ncol=1,
                   align = "v",
                   rel_heights = c(1, 1.3))

ggsave(plot = Figure_4, "outputs/Figure_4.png", dpi = 600, height =40, width = 30, units = "cm")

